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Abstract 

A general framework for the simulation of reaction-diffusion sys- 
tems with probabilistic cellular automata is presented. The basic re- 
action probabilities of the chemical model translate directly into the 
transition rules of the automaton, thus allowing a clear comparison 
between simulation results and analytic calculations. This framework 
is then applied to simulations of hypercycle-systems in up to three di- 
mensions. 

Furthermore, a new measurement quantity is introduced and applied to 
the hypercycle-systems in two and three dimensions. It can be shown 
that this quantity can be interpreted as a measure for the macroscopic 
order of the hypercycle systems. 



1 The automaton 



In the past a lot of work was done in the field of algorithmic chemistry |T|, ^ 
and chemistry on cellular automata (CA) Our goal is to give a com- 

plete framework which begins with the determination of the transition rules 
from basic reaction probabilities and ends with the comparison between the 
emergence of spatially ordered structures and a macroscopic measurement 
quantity in hypercycle-systems. 

In this paper we can only give a brief introduction of our probabilistic CA 
with asynchronous updating. Our aim was to establish a one-to-one corre- 
spondence between real chemical reactions (defined by their reaction prob- 
abilities) and the simulations on a CA (defined by the transition rules). 
Therefore, we explain the calculation of these transition rules for a simple 
particle production process in detail and then present the transition rules 
for more complex reaction types without derivation. 
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1.1 Particle production processes 

In our system, particles can be produced through spontaneous or induced 
reactions. The result of both reactions is the change of an empty cell to an 
occupied one. Since we evaluate production processes for empty cells, we 
will usually have the situation that several different reactions are possible, 
resulting in a different final state for the previously empty cell. If we think, 
e.g. of spontaneous self-replication, all the nearest neighbors of the empty 
cell are in principle equally likely to reproduce into this cell. In chemistry 
such symmetry breaking phenomena are induced by statistical fluctuations, 
e.g. temperature gradients. 

Let Pi, i A, he the probability for a single species i to replicate into the 
empty cell by a special reaction, e.g. self-replication. Here A denotes an 
index set containing all indices which are necessary to enumerate all species 
in the given context and Ni will denote the number of molecules of species 
i. 

In the case there are several nearest neighbor cells occupied, we first have 
to calculate the total probability that at least one arbitrary nearest- 
neighbor-molecule will occupy the empty cell. We get 

p*°*=i-n(i-K)'^^ (1) 

From we can calculate the effective probabilities for the nearest neighbor 
molecules to replicate in the given arrangement of species. We will denote 
these probabilities with p^^. In this calculation we have to make sure that 
the ratio Pi/pj is equal to pf^/pf^, G A. Since the sum over all pi is 
not necessarily 1, we introduce a normalization factor X^ieyt-^iK-Q We get 

„tot 

pf=p,^^. (2) 



Note that is the sum of Nipl^ over all i. 

1.1.1 Type I reactions 

Reactions of type I, 



represent the group of spontaneous reactions with two reaction products. 
After some work and with the calculations done above we get: 

„eff _ „ 1 ~ ritg^Cl ~ ■P^,r?^,n)^' 

Pi,m,n ~ Pi,m,n \^ AT n ' ^ ' 



^The normalization factor is introduced to ensure that the sum over aU p"*^ and 1 - 
is 1. 
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The corresponding transition rule then reads: 

A cell with state zero will change to state n and a nearest-neighbor cell with 



state i will change to m with the probability p':^ 



t,m,n' 



1.1.2 Type II reactions 

Reactions of type II, 



describe a group of two-body reactions with three reaction products. Obvi- 
ously, such reaction processes are more difficult to handle since we have to 
implement the determination of the correct reaction partner. The related 
quantities are primed. 
After some calculations we get: 



eflf ^ / 1 ntg^(l PiJ,m,n,o) ' / 

Pi,j,m,n,o Pi,j,m,n,o at / • 



with 



m 

P'i,3,m,n,o = P{^' ^ 0) j;^Pi,j,m,n,o- (5) 

p{N' 7^ 0) denotes the probability that there is at least one molecule to 
react with. For the derivation of Pij^m,n,o have implicitly used the idea 
of an attractive force between the molecules which means that the reactive 
molecule therefore "jumps" to a reaction partner. On the other hand, if we 
do not assume such an attractive force, which means a molecule can also 
"jump" into a hole, we have to substitute Pij^m,n,o with iVj/-^max- 
Then we can write down the transition rule: 

The probability that a cell in state zero, with a nearest-neighbor cell in state 

i and with a "reaction" cell in state j will change to o, while the nearest- 
neighbor cell changes to m and the "reaction" cell changes to n, is p^^ 



i,j,m,n,o- 



1.1.3 Particle decay 

The decay process. 



can change an occupied cell into an empty one. The transition rule descrip- 
tion is rather easy: 

The probability that a cell in state i will change to a cell in state zero is pf^'^ . 
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1.2 Diffusion 



The diffusion is an additional process. Since we use an asynclironous au- 
tomaton, an often used diffusion algoritlim for CA |5[ is not necessary. 
Instead, we use a more natural algorithm that can be divided into three 
subprocesses: 

1. An arbitrary cell of the CA is chosen. 

2. One of it's nearest neighbor cells is determined at random. 

3. The states of the two cells are interchanged. 

1.3 Automaton update 

We use asynchronous updating, since this seems to be the more natural way 
1^. The basic processes are then executed in the following way: 
Choose an arbitrary cell at random. If the cell is occupied we use the decay- 
transition rule. Otherwise the transition rule for reactions of type I or of 
type II or a combination of both are used. 

The reaction processes and the diffusion process are then combined to the 
full simulation in the following way. Let k be the total number of cells in 
a CA and d the ratio of the number of basic processes and the number of 
diffusion steps. Then a CA-time-step t consists of k basic processes and dk 
diffusion steps. It is important that the basic processes and the diffusion 
steps are mixed up very well. 

Note that t is independent of the system's size, since we scale the number 
internal processes with it. 



2 Simulation 

We have simulated systems with (symmetrical) hypercycles |^, ^, which 
means we set 

^^-'^ = jo else 

and 



jpcat j£ ^ _ _ J- _ ^ g^y^^ ^ _ J- _j_ £ (cyclic), 

Pi,j,m,n,o ~ \ Q Y 
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2.1 Stable patterns 

In this part of the paper we will show the well-known spiral waves Q in 
two dimensions and the emergence of a stable scroll wave structures in 
three dimensions. In the latter case we have used a two-dimensional spiral 
as initial condition. We will see the system going through different unstable 
patterns. For reasons of clear visualization, in the pictures we show only 
one of the six species present in the simulation. 

The three-dimensional analog of a spiral wave is a scroll wave Q. As we 
can see in Fig. ^, even in the two dimensional case, structures can be very 
complex. Since there are more degrees of freedom in three spatial dimen- 
sions, the structures arising from non-ordered initial conditions are much 
more complex. Therefore, we use a 2d spiral as a symmetry-breaking initial 
condition for our simulations to simplify the emergent structures. 
We can see the spiral's growth in the z-direction (see Fig. ||, ^). As a matter 
of fact, the reproduction activity of a species becomes high if it is adjacent to 
it's "supporter". Therefore, in a layer above, molecules are first created over 
the small region where their own and their "supporter" are present. Since 
the spirals rotate in the "supporter" direction, we can see the phenomenon 
that the new spirals, growing on top of the old ones, are phase shifted in 
the rotation direction. This continues until the space in z-direction is oc- 
cupied. The resulting structure looks like a rotating screw. When all the 
available space is occupied, the molecules cannot replicate freely any more. 
Since space is rare than, the molecules have to compete for empty space to 
replicate. For this reason, the molecules begin to arrange themselves in a 
different way so that the "free" screw begins to change to a cone structure. 
This cone structure is stable so that a real scroll wave cannot be seen. 
The whole re- arrangement process takes much more time than the occupa- 
tion at the beginning. The pictures show the emergence of the cone struc- 
ture. At the top of the cone the emergence of other small cones can be seen 
so that the whole structure becomes more complex. 

2.2 A microscopic measurement quantity 

Now we will introduce a measurement quantity with which we can measure 
the arrangement of the molecules in our automaton. The definition of this 
order-quantity, called a^, is rather easy: 

For every molecule, we take the nearest neighbor molecules. If an adjacent 
molecule is a "supporter", the variable S is incremented by one. Otherwise 
the variable is incremented by one. If a cell is not occupied, nothing is 
done. With these microscopic values we can define the macroscopic quantity 
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Figure 1: The first two pictures show the system's behavior shortly after the 
"injection" of the 2d spiral. We can see the emergence of a rotating screw. 
The left picture was taken at t = 150 and the right one at t = 250. 




Figure 2: If the space is nearly completely occupied the molecules cannot 
replicate freely and have to compete for empty space. Therefore the struc- 
ture changes from a screw to a cone-like pattern. 
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Figure 3: The three pictures show the time dependent behavior of the sys- 
tem. The times when the pictures have been taken are (from left to right): 
t = 0,t = 120 and t = 9000. See also Fig. ^ and text for further information. 



= WTks' 

where /c > is a tuning constant. The larger k, the more sensitive is on 
5. 

Since both N,S >0, the quantity ranges between and 1. We will later 
see that Ok is even constant for large t. (Obviously, if is constant, Uk', 
with another k' ^ k, is also constant.) 

Moreover, since we do not take into account empty cells, is independent 
of the system's concentration and therefore independent ofp'^'^'^. 
For our simulations we have used k = 10 and applied ctio to both the two 
and three dimensional case. 

Note that there are a lot of possible microscopic/macroscopic quantities 
which all can describe the system from different points of view. 



2.3 Two dimensions 

For the two dimensional case we have used two different initial conditions. 
In both cases the parameter set is p'^'^'^ = 0.1, p^'^^^ = 0.05, p^^* = 0.9 and 
d = l. 



First, the hypercycle-system was initialized at random, see Fig. |. We 
have seen aiQ converging rapidly to a constant value at the very beginning. 
After t = 120, however, aio changes only in very small amounts. In compar- 
ison, the global structure of the system also changes very slowly, see Fig. ^. 
It seems that the system's rough structure is determined at the very begin- 
ning and will then only be refined, later on. This is very interesting, since 
there are a lot of possible stable patterns. However, the decision is made 
very early. 

For the asymptotic value of aio we have found aio ~ 0.72, independent 
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Figure 4: The picture shows the concentrations of the six hypercycle species 
(soHd Unes) and the value of aio versus time for a system with random initial 
condition. 




Figure 5: The three pictures show the time dependent behavior of the sys- 
tem. The times when the pictures have been taken are (from left to right): 
t = 0, t = 140 and t = 550. See also Fig. ^ and text for further information. 

of the number of species in the hypercycle and of the system's size. Fur- 
thermore, ctfc depends only weakly on d for d < 1. The next step was to 
examine a system where only one stable final pattern exits. Therefore, we 
have used a special initial condition which is displayed in Fig. |5[ The final 
structure will be of course a single rotating spiral. 

What we have seen was aio changing comparably to the transformation of 
the triangular structure to the spiral structure. Since the major part of this 
transformation is done at the beginning, a also increases fast during the first 
simulation steps and then increases slower and slower. The system behaves 
similar to the case with random initialization, although the convergence of 
ak to it's asymptotic value is significantly slower. But what is more surpris- 
ing, the asymptotic value of qiq is again the same, aio ~ 0.72. 
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Figure 6: The picture shows the concentrations of the six hypercycle species 
(sohd Unes) and the value of aio versus time for a system with ordered 
triangles as initial condition. 



We will now go a step beyond and apply aio on our three dimensional sim- 
ulator. 

2.4 Three dimensions 

As we have seen above (Fig. [l|, ^ in three dimensions there are several 

different unstable structures, beginning with a single two dimensional spiral 

as initial condition and ending with a stable cone-structure. 

In Fig. 1^ we can see aio starting with aio ~ 0.72. But then we see a^ 

increasing in the time interval when the spiral occupies the empty space by 

growing in the z-direction. After that, when the available space is occupied, 

new boundaries are present and aio decreases to the former value. 

This is a bit surprising since the cone structures are quite different and at the 

first glance there seems to be no reason why the macroscopic order should 

be similar. 



3 Conclusions 

In this paper we have shown pictures of the transformation from two di- 
mensional patterns to three dimensional ones and have therefore seen the 
system walking through different states of arrangement. But what is more 
interesting is the fact that our macroscopic quantity a^ (here: aio) is quite 
constant after a short time and has even the same asymptotic value in dif- 
ferent settings. 

In all our simulations, especially in two dimensions, we have seen aio con- 
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Figure 7: The picture shows the concentrations of the six hypercycle species 
(sohd Unes) and the value of aio versus time for a system in three spatial 
dimensions. Note that a single two dimensional spiral was used as initial 
condition. 



verging rapidly to the fixed value, ckiq ~ 0.72. Thereof we can deduce that 
such catalytic systems have the urge to determine their final structures at 
the very beginning. Later on structures are only refined. 
Although the hypercycles have been studied very well in the past years some 
new features can still be found. We think that the definition of microscopic 
quantities and their interpretation in a macroscopic sense is one of the right 
ways to learn more about catalytic systems which can be even more complex 
than such "simple" symmetric hypercycles. 



4 Outlook 

It seems to be interesting for future work to study whether stability against 
parasites is possible in three dimensions, similar to the two-dimensional case 
IC]. Furthermore, the behavior of Uk in such parasite-systems might be 
interesting. 
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